% This code tests for structual break.
close all; clear; clc;

fig_save = 'figE1_chow_test.png';

%% load saved output from 'est1' and 'est2'
% If the .mat file is not in folder, run 'est1_estimation' with Tend = 2015
load est1_predict_HSSB_2015.mat
year = 1950:2015;

% If the .mat file is not in folder, run 'est2_regression1' with Tend = 2015
load est2_output1_residual_2015.mat

% If the .mat file is not in folder, run 'est2_regression2' with t1 = 1996
load est2_output2_residual_1996.mat

%% compute residual of H and SSB (unit: kg)
residH = (H - H_hat);
residSSB = (SSB2-SSB_hat);

%% Chow test -- residual of H
bp = 6:length(year)-5; yearTest1 = year(bp);
[h1, pValue1,stat1,cValue1] = chowtest(residH(1:end-1),residH(2:end),bp,...
    'Display','off');
fprintf('\n\nTest break(s) in H residuals:\n')
disp(yearTest1(h1))

%% Chow test -- residual of SSB
residSSB2 = residSSB(~isnan(residSSB));
yearTest2 = year(~isnan(residSSB));
bp = 6:length(yearTest2)-5; yearTest2 = yearTest2(bp);
[h2, pValue2,stat2,cValue2] = chowtest(residSSB2(1:end-1),residSSB2(2:end),bp,...
    'Display','off');
fprintf('\n\nTest break(s) in SSB residuals:\n')
disp(yearTest2(h2))

%% Chow test -- residual of regression (constant cost)
bp = 6:length(year)-5; yearTest3 = year(bp);
[h3, pValue3,stat3,cValue3] = chowtest(resid1(1:end-1),resid1(2:end),bp,...
    'Display','off');
fprintf('\n\nTest break(s) in Regression1 residuals:\n')
disp(yearTest3(h3))

%% Chow test -- residual of regression (constant cost)
bp = 6:length(year)-5; yearTest4 = year(bp);
[h4, pValue4,stat4,cValue4] = chowtest(resid2(1:end-1),resid2(2:end),bp,...
    'Display','off');
fprintf('\n\nTest break(s) in Regression2 residuals:\n')
disp(yearTest4(h4))

%% Plot pvalue of Chow test
tlt1 = sprintf('(a) For residuals of H(t),\n null hypothesis of Chow test cannot be rejected');
tlt2 = sprintf('(b) For residuals of SSB(t),\n null hypothesis of Chow test cannot be rejected');
tlt3 = sprintf('(c) For residuals of OLS (constant cost),\n null hypothesis of Chow test cannot be rejected');
tlt4 = sprintf('(d) For residuals of OLS (cost shift in 1996),\n null hypothesis of Chow test cannot be rejected');

fig3 = figure('Position',[0 0 900 500]);
subplot(2,2,1) % -- H --
hold on
plot(yearTest1,pValue1,'LineWidth',1.5)
scatter(yearTest1(h1),pValue1(h1),'r*')
plot(year,0.05+0*year,'k:','LineWidth',2)
ylim([0 1])
xlim([year(46) year(end-5)])
title(tlt1) %'Chow test on residuals of H(t)'
xlabel 'year'
ylabel 'p-value'

subplot(2,2,2) % -- SSB --
hold on
plot(yearTest2,pValue2,'LineWidth',1.5)
scatter(yearTest2(h2),pValue2(h2),'r*')
plot(year,0.05+0*year,'k:','LineWidth',2)
ylim([0 1])
xlim([year(46) year(end-5)])
title(tlt2) %'Chow test on residuals of SSB(t)'
xlabel 'year'
ylabel 'p-value'

subplot(2,2,3) % -- reg1 --
hold on
plot(yearTest3,pValue3,'LineWidth',1.5)
scatter(yearTest3(h3),pValue3(h3),'r*')
plot(year,0.05+0*year,'k:','LineWidth',2)
ylim([0 1])
xlim([year(46) year(end-5)])
title(tlt3) %'Chow test on residuals of OLS regression with constant cost'
xlabel 'year'
ylabel 'p-value'

subplot(2,2,4) % -- reg2 --
hold on
plot(yearTest4,pValue4,'LineWidth',1.5)
scatter(yearTest4(h4),pValue4(h4),'r*')
plot(year,0.05+0*year,'k:','LineWidth',2)
ylim([0 1])
xlim([year(46) year(end-5)])
title(tlt4) %'Chow test on residuals of OLS regression with cost shift in 1996'
xlabel 'year'
ylabel 'p-value'


saveas(fig3,fig_save)